#define M_PI       3.14159265358979323846

#include <vector>
#include <complex>
#include <cmath>
#include <iostream>

using Complex = std::complex<float>;
using IQVector = std::vector<Complex>;

class IQResampler {
public:
    // 设计线性相位低通滤波器
    static std::vector<float> design_lpf(int taps, float cutoff) {
        std::vector<float> coeffs(taps);
        float sum = 0.0f;

        // Sinc函数加汉宁窗
        for (int n = 0; n < taps; ++n) {
            float x = n - (taps - 1) / 2.0f;
            coeffs[n] = (x == 0) ? 2.0f * cutoff
                : std::sin(2 * M_PI * cutoff * x) / (M_PI * x);
            coeffs[n] *= 0.5f - 0.5f * std::cos(2 * M_PI * n / (taps - 1));
            sum += coeffs[n];
        }

        // 归一化保证增益=上采样倍数
        for (auto& c : coeffs) c *= 1.0f / sum;
        return coeffs;
    }

    // 修正后的上采样（功率精确保持）
    static IQVector upsample(const IQVector& input, int factor) {
        if (factor <= 1) return input;

        // 零填充
        IQVector output(input.size() * factor, Complex(0, 0));
        for (size_t i = 0; i < input.size(); ++i) {
            output[i * factor] = input[i];
        }

        // 设计滤波器（增益补偿factor倍）
        auto lpf = design_lpf(64, 0.5f / factor);

        // 应用滤波器（含功率补偿）
        apply_filter(output, lpf, factor);
        return output;
    }

    // 下采样（保持原实现）
    static IQVector downsample(const IQVector& input, int factor) {
        if (factor <= 1) return input;
        auto lpf = design_lpf(128, 0.5f / factor);
        IQVector filtered = input;
        apply_filter(filtered, lpf, 1.0f);

        IQVector output;
        for (size_t i = 0; i < filtered.size(); i += factor) {
            output.push_back(filtered[i]);
        }
        return output;
    }

private:
    // 改进的滤波器应用（自动补偿功率）
    static void apply_filter(IQVector& signal,
        const std::vector<float>& coeffs,
        float compensation_gain = 1.0f) 
    {
        const int half = coeffs.size() / 2;
        IQVector temp(signal.size());

        for (size_t i = half; i < signal.size() - half; ++i) {
            Complex sum(0, 0);
            for (size_t j = 0; j < coeffs.size(); ++j) {
                sum += signal[i + j - half] * coeffs[j];
            }
            temp[i] = sum * compensation_gain;
        }

        // 边界零填充
        for (size_t i = 0; i < half; ++i) {
            temp[i] = Complex(0, 0);
            temp[signal.size() - 1 - i] = Complex(0, 0);
        }

        signal = temp;
    }
};

// 功率测试工具
void test_power() {
    const int N = 10000;
    IQVector iq_data(N, Complex(1.0f, 1.0f)); // 理论功率=2.0

    std::cout << "=== 功率补偿验证 ===\n";
    std::cout << "输入功率: " << 10 * log10(2.0) << " dB\n";

    for (int factor : {2, 4, 8}) {
        auto upsampled = IQResampler::upsample(iq_data, factor);
        float power = 0;
        for (auto& s : upsampled) power += norm(s);
        power /= upsampled.size();
        std::cout << "上采样 " << factor << "x 功率: "
            << 10 * log10(power) << " dB (误差: "
            << 10 * log10(power / 2.0f) << " dB)\n";
    }

    for (int factor : {2, 3, 4}) {
        auto downsampled = IQResampler::downsample(iq_data, factor);
        float power = 0;
        for (auto& s : downsampled) power += norm(s);
        power /= downsampled.size();
        std::cout << "下采样 " << factor << "x 功率: "
            << 10 * log10(power) << " dB (误差: "
            << 10 * log10(power / 2.0f) << " dB)\n";
    }
}

int test_resample_main() {
    // 运行功率测试
    test_power();

    // 示例使用
    IQVector input(100, Complex(1.0f, 0.5f));
    auto upsampled = IQResampler::upsample(input, 2);
    auto downsampled = IQResampler::downsample(input, 3);

    std::cout << "\n=== 运行示例 ===\n";
    std::cout << "上采样后长度: " << upsampled.size() << "\n";
    std::cout << "下采样后长度: " << downsampled.size() << "\n";

    return 0;
}